#------火山图-----
load("ad_ct.RData")
table(pro_data$group)
add_file=add_file[-35,]
data1=cbind(add_file,pro_data)



#https://www.jianshu.com/p/e2828ef9c7e5
#加载包
library(ggplot2)
library(ggrepel)
#abs 绝对值函数
foldchange=2^ad_ct$log2FC    #log2的反运算
#分三种情况填入新列
ad_ct$threshold = (ifelse(ad_ct$P < 0.05 & foldchange > 1.5, 'UP',
                          ifelse(ad_ct$P < 0.05 &foldchange<= 1/1.5 ,
                                 'Down','NoSignifi')))
#采用ifelse的嵌套添加新列，标识基因是否特异表达(刚学的，觉得挺牛)

ad_ct$ID2=(ifelse(ad_ct$threshold != "NoSignifi",ad_ct$ID,''))
#添加ID2列，单独记录上行表达和下行表达的基因名，


table(ad_ct$threshold)  
#table查看发现是down、nosignifi、up的顺序，指示后来颜色的命名
ggplot(ad_ct,aes(x=log2FC,y=-log10(P),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("green","gray","orange"))+#确定点的颜色
  xlab('log2 (FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-log2(1.5),log2(1.5)),lty=3,col="black",lwd=1) +
  #添加竖线|FoldChange|>1.5
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=1)+
  #添加横线P<0.05 
  geom_text_repel(aes(label=ad_ct$ID2),color="black",size=3)
#添加基因标签

#ADdata2[==0]=NA  #去0,将0转换成NA





#②
#abs 绝对值函数
load("as_ct.RData")
foldchange=2^as_ct$log2FC    #log2的反运算
#分三种情况填入新列
as_ct$threshold = (ifelse(as_ct$P < 0.05 & foldchange >=1.5, 'UP',
                          ifelse(as_ct$P < 0.05 &foldchange<= 1/1.5 ,
                                 'Down','NoSignifi')))
#采用ifelse的嵌套添加新列，标识基因是否特异表达(刚学的，觉得挺牛)

as_ct$ID2=(ifelse(as_ct$threshold != "NoSignifi",as_ct$ID,''))
#添加ID2列，单独记录上行表达和下行表达的基因名，


table(as_ct$threshold)  
#table查看发现是down、nosignifi、up的顺序，指示后来颜色的命名
ggplot(as_ct,aes(x=log2FC,y=-log10(P),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("green","gray","orange"))+#确定点的颜色
  xlab('log2 (FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-log2(1.5),log2(1.5)),lty=3,col="black",lwd=1) +
  #添加竖线|FoldChange|>1.5
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=1)+
  #添加横线P<0.05 
  geom_text_repel(aes(label=as_ct$ID2),color="black",size=3)

#添加基因标签

#ADdata2[==0]=NA  #去0,将0转换成NA


  #③
  #abs 绝对值函数
load("mci_ct.RData")
foldchange=2^mci_ct$log2FC    #log2的反运算
#分三种情况填入新列
mci_ct$threshold = (ifelse(mci_ct$P < 0.05 & foldchange >=1.5, 'UP',
                          ifelse(mci_ct$P < 0.05 &foldchange<= 1/1.5 ,
                                 'Down','NoSignifi')))
#采用ifelse的嵌套添加新列，标识基因是否特异表达(刚学的，觉得挺牛)

mci_ct$ID2=(ifelse(mci_ct$threshold != "NoSignifi",mci_ct$ID,''))
#添加ID2列，单独记录上行表达和下行表达的基因名，


table(mci_ct$threshold)  
#table查看发现是down、nosignifi、up的顺序，指示后来颜色的命名
ggplot(mci_ct,aes(x=log2FC,y=-log10(P),color=threshold))+
  geom_point()+
  scale_color_manual(values=c("green","gray","orange"))+#确定点的颜色
  xlab('log2 (FoldChange)')+#修改x轴名称
  geom_vline(xintercept=c(-log2(1.5),log2(1.5)),lty=3,col="black",lwd=1) +
  #添加竖线|FoldChange|>1.5
  geom_hline(yintercept = -log10(0.05),lty=3,col="black",lwd=1)+
  #添加横线P<0.05 
  geom_text_repel(aes(label=mci_ct$ID2),color="black",size=3)
  
  #添加基因标签





#--------提取筛选后的相关基因列表------
label=ad_ct$threshold
up=label!="NoSignifi"
ad_change_protein <- ad_ct[up,]
save(ad_change_protein,file = "ad_change_protein.RData")

#-——-整合出差异蛋白信息----
load("as_change_protein.RData")    #as组差异表达蛋白
load("ad_change_protein.RData")    #ad组差异表达蛋白
load("mci_change_protein.RData")   #mci组差异表达蛋白
names(ad_change_protein)[4]=names(as_change_protein)[4]#统一列名
ad_change_protein=ad_change_protein[,-5]    #删除多余列
change_protein1=rbind(ad_change_protein,as_change_protein,mci_change_protein)#整合
change_protein=unique(change_protein1$ID)    #共得到115个候选蛋白


load("BLSA_data.RData")          #调取转置前总文件      
pos <- blsa_pro$Uniprot.ID %in% change_protein
blsa_pro <- blsa_pro[pos,]
rownames(blsa_pro)=blsa_pro[,3]
blsa_pro.t2 <- as.data.frame(t(blsa_pro[,-1:-3])) #删除冗余列，转置
#整合完毕

#--------树状图-----
#  树状图
dist_ADdata=dist(blsa_pro.t2,method="manhattan")  
#曼哈顿法求距离，每一行之间的距离
dist_ADdata=hclust(dist_ADdata,method="average")
#hculst,根据距离分类,有average,centroid,median多种
plot(dist_ADdata,hang=-1,cex=.8,col=c("red","lightblue","green","orange"))  
#hang=-1会对齐最下方  ，cex.8是点的大小

install.packages("dendextend")
library(dendextend)
class(dist_ADdata)
#聚类结果原始“hclust”类先转化为"dendrogram"
dend <- as.dendrogram(dist_ADdata)
#dendextend程序包的函数对不同的分组进行标色
dend %>% set("branches_k_color") %>% plot

# 使用标准色显示聚类簇
clusters <- cutree(dend,5)[order.dendrogram(dend)]
dend %>% set("branches_k_color", k=5,
            value = unique(clusters) + 1) %>% plot
# 增加颜色条带
colored_bars(clusters + 1,
             y_shift = -0.5,
             rowLabels = paste(5, "clusters"))



#---------直方图----------
library(psych)
add_file.des <- describe(add_file,IQR=TRUE)
hist(unlist(blsa_pro.t2),
     col=rainbow(15),breaks=c(min(blsa_pro.t2),-2,-1,-0.5,-0.2,0,0.2,0.5,1,2,
                              max(blsa_pro[,-1:-3])),main="数据优化后的直方图",
     xlab = "MFG EXPRESSION")
lines(density(unlist(blsa_pro.t2))) #密度曲线

#-----------基本特征描绘——————

library(psych)
add_file.des <- describe(add_file)
#blsa_pro:MFG LOG2 EXPRESSION MATRIX
blsa_pro=blsa_pro[,-1:-3]
blsa_pro=abs(blsa_pro)             #取绝对值，防止效应抵消
blsa_pro.des <- describe(blsa_pro)   #得到针对表达量的数据基本特征

a=a[,-4:-5]                   #这几步为只保留基因名和均值列
a[,2]=row.names(a)
a[,1]=gsub("_[0-9]+","",row.names(a))
library(openxlsx)
write.xlsx(a,"outcome1.xlsx")

oneway.test(a[,3]~a[,1])       #四种分型间的分析   
#F = 3.3342, num df = 3.000, denom df = 27.295, 
#p-value =0.03399
#说明有显著差异
a[,1][a[,1]=="AS"]="AD"      #比较AD和AS的整合组与ct的区别
b=a[-48:-58,]                  #48-58为mci组数据，删除 
oneway.test(b[,3]~b[,1])       #两种分型间的分析   
#F = 7.5104, num df = 1.000, denom df = 23.516, 
#p-value =0.01152
#说明有显著差异










